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FOREWORD 


The  principal  purpose  of  this  paper  is  to  present  methods  for  determinir.g 
minimum-time  ship  routes  on  a  digital  computer.  The  paper  was  written  in 
two  parts.  The  first  treats  in  detail  a  numerical  method  for  determining  the 
course  which  requires  minimum  time  to  go  from  one  specified  place  to  another 
when  the  speed  of  the  ship  is  a  known  function  of  position,  heading,  and  time. 

The  second  part  of  the  paper  treats  various  related  topics.  The  first  is 
the  corresponding  problem  in  the  simpler  case  wherein  the  speed  is  a  functic'n 
of  position  and  heading  but  not  of  time.  The  second  is  a  discussion  of  metho<is 
for  obtaining  the  curves,  called  equal-time  curves  or  isochrones,  which  yield 
the  maximum  distance  the  ship  can  attain  at  any  particular  time  by  choosing 
various  courses.  These  are  discussed  for  both  the  stationary  and  time- 
varying  speed  fields;  they  are  rather  important  because  they  are  easily 
understood  and  interpreted,  and  closely  parallel  numerical  methods  now  in  use. 
Third,  the  problem  of  effecting  rendezvous  between  two  ships  is  treated. 
Finally  a  brief  discussion  is  given  of  various  problems  and  difficulaties  which 
may  be  encountered  in  computation. 

It  has  been  the  intention  particularly  in  the  first  part  to  give  the  numerical 
routines  and  associated  discussions  in  sufficient  detail  that  a  person  familiar 
with  numerical  methods  and  computers  could  immediately  program  and  run 
the  problem.  The  method  is  based  on  procedures  which  G.  A.  Bliss  introduced 
in  ballistics  for  calculating  differentials.  These  are  applied  in  a  Newton- 
Raphson  iteration  to  determine  the  course.  It  is  the  author's  opinion  that  the 
major  outstanding  part  of  the  ship- routing  problem  is  the  collection  of  re¬ 
liable  empirical  data  describing  the  speed  of  the  ship. 

It  is  hoped  that  this  paper  will  serve  as  an  introduction  to  the  adjoint  sys¬ 
tem,  calculus  of  variations,  and  optimum  control  theory  along  lines  which  aie 
currently  being  actively  pursued  in  this  country  and  in  Russia,  including  a 
method  of  solution.  The  first  part  of  this  paper  has  been  accepted  for  publica¬ 
tion  in  NAVIGATION,  Journal  of  the  Institute  of  Navigation. 
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A  GENERAL  NUMERICAL  METHOD  FOR  DETERMINING  OPTIMUM  SHIP 

ROUTES 

Frank  Do  Faulkner* 

A  method  is  given  for  determining  optimum^  ship  routes  on  a 
digital  computero  This  paper  is  limited  to  the  problem  of  deter- 
mlni^mlnlmum  tlme„cour,ses^when  the  speed  is  a  known  function  of 
the  position,  heading,  and  time. 

The  method  is  based  on  Bliss's  methods  for  calculating  dif¬ 
ferentials,  usually  applied  In  a  Newton  Iteration  to  determine 
the  course.  It  Is  very  general,  but  It  requires  a  knowledge  of 
Bliss's  adjoint  methods.  The  paper  Is  selfoontalned;  a  simple 
case  Is  given  In  detail  and  others  outlined.  The  advent  of  high¬ 
speed  digital  computers  In  the  past  few  years  now  makes  the  solu- 
tion  of  such  problems  feasible. 


1.  Statement  of  problem. 

shall  take  the  equations  of  motion  to  be  expressed  in  the 

form 

o 

=  V  cos  p 
ly  =  V  sin  p 

where  x,y  are  position  coordinates,  p  Is  a  control  variable, 

''  =  ''(x,y,p,t)  Is  a  known  function  of  period  211  in  p  with 
continuous  derivatives,  t  Is  time,  and  the  dot  C)  over  a 
variable  Indicates  Its  time  derivative.  If  x,y  are  rectang¬ 
ular  coordinates,  then  v  Is  the  speed  and  p  la  the  angle 
between  the  velocity  vector  and  the  x  axis  ,  the  heading  angle; 
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If  they  are,  say,  latitude  and  longitude,  then  there  is  no  such 
simple  interpretation. 

We  will  consider  first  the  problem  of  going  from  one  specified 
point,  (0,0),  to  another,  (X,Y) ,  with  minimum  time  T.  The  problem 
is  to  choose  the  control  variable  p  as  a  function  of  time  to  ef¬ 
fect  this,  and  to  determine  the  corresponding  curve  in  x,y,t  space. 
In  a  later  section  general  conditions  will  be  given  for  minimum 
time  courses,  as  in  rendezvous,  etc.,  and  corner  conditions  which 
would  be  of  interest  for  a  sailboat  which  must  tack. 


2.  Variational  equations,  Euler  equations. 

In  this  section  some  formulas  for  differentials  are  derived, 

by  procedures  which  G.  A.  Bliss  introduced  in  Ballistics.  Let  us 

consider  two  neighboring  paths,  whereon  the  values  of  p  differ 

r  T,  , 


by  an  amount  6p;  it  is  assumed  that 


6p|  dt  is  small, 


Then 


the  change  5x,6y  in  x,y  satisfy  the  variational  equations 


(2) 


'6x  =  cos  p  +  Vy6y  cos  p  +  (VpCOS  p  -  v  sin  p)6p 

l6y  =  sin  p  +  v^by  sin  p  +  (v^sin  p  +  v  cos  p)6p; 


subscripts  here  Indicate  partial  derivatives,  v  =  dv/dx,  etc. 

Let  us  multiply  these  two  equations  through  by  two  new 
variables  (unspecified  at  present  but  to  be  identified  a^ 

some  stage  as  Lagrange  multipliers),  collect  terms,  and  Integrate 
to  get 
(3) 


rT 


[)v(6x  -  V  6x  cos  p  -  V  by  cos  p  -  V  6p  cos  p  +  v6p  sin  p) 

0  X  y  p 

h  ^(by  -  v^6x  sin  p  -  Vy.6y  sin  p  -  v^bp  sin  p  -  vbp  cos  p)]d 


=  0 

We  may  integrate  this  by  parts  and  rewrite  it  as 
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r  T 

(4)  [9i6x  +y-'£'y]o  =  [6x(^+Xv^cos  p  +y^>'v^sln  p)  + 

+  &y(y^  +  XVyCOS  p  +yJVySln  p) 

+  6p(A{v^  ^cos  p  -  V  sin  p}-  +yO'{VpSln  p  +  v  cos  p}')]dt 

ti+ 


-  [v(Xcos  p  +y^sin  p)]  6t^ 


where  is  a  symbol  for  any  point^.or  points  where  p  is  dis¬ 

continuous  as  a  function  of  time,  a  steering  corner.  Probably 
there  will  be  no  such  point  for  ships  and  we  shall  drop  the  last 
term  for  the  present.  To  simplify  this  equation,  let  us  choose 
as  solutions  to  the  system 


(5) 


■X  +>  V  cos  p  +  /-^v  sin  p 

X  X 

* 

yu*  +  XVyCOS  P  +  ylv'VySin  p 


0 

0  , 


so  that  the  coefficients  of  5x,6y  in  (4)  vanish.  This  defines 
the  system  (5)  which  is  adjoint  to  the  system  (3)o 

In  the  case  of  most  Interest,  the  initial  values  of  x,y  are 
assumed  known;  then  6x(0)  =  0  =  6y(0)  ,  and  (4)  becomes 


(6)  0y(T)6x(T)  +  y(^(T)6y(T)  = 

cos  p  -  V  sin  p)  +y»-'(VpSln  p  +  v  cos  p)]6p  dt 


TT- 

=  A^v  6p  dt, 

Jo  P 

-A.  ^ 

where  A  =  "Xi  +^3*  and  "v  =  v(l  cos  p  +  j  sin  p).  This  is  the 
fundamental  differential  formula  connecting  the  change  in  end 
values  of  x,y  with  the  variation  of  the  control  variable. 

It  is  a  very  general  condition  for  ’any  optimum  that  there  is 
some  solution  A  to  the  adjoint  system,  to  be  found  in  the  solu¬ 
tion  of  the  problem,  and  for  that  solution,  ^  must  satisfy  the 
condition 
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(7) 


A’v  =  extreme, 


as  a  function  of  the  control  variable  p  at  every  point  of  the 
course.  This  condition  implies  that  the  coefficient  of  6p  in 
(6)  vanish 


for  that  solution.  Comment.  The  proof  of  this  condition  is  es¬ 
sentially  the  proof  of  the  fundamental  lemma  in  the  calculus  of 
variations  (see  Courant*[l]  p  200).  Cases  wherein  it  does  not 
determine  a  unique  path  are  easily  constructed,  but  seem  to  be  of 
little  practical  significance  and  will  not  be  considered  further. 

Equations  (5) »  (8)  are  called  the  Euler  equations  in  calculus 
of  variations.  The  curves  whereon  equations  (1),(5)»(8)  are  satis 
fled  are  called  extremals  of  the  family  (1).  Equation  (7)  or 
(8)  may  be  replaced  by  the  condition  that  p  be  chosen  so  that 


(9) 


0 


this  approach  has  received  considerable  attention  recently  through 
the  work  of  the  Russian  mathematician  Pontrlagin  [2]  on  control. 

The  various  philosophies  of  approach  represented  by  (7) » 

(8), (9)  are  equivalent  for  this  problem,  and  none  eliminates 
the  basic  problem  of  solution;  namely,  conditions  are  given  at 
two  values  of  t,  (0,T),  with  the  second  one  unknown,  and  a  con¬ 
stant  is  needed  to  integrate  equations  (1),(5)»(8).  In  the  next 
sections  a  method  of  solution  is  given. 

3.  Differential  formulas. 

Let  us  consider  a  fundamental  set  of  solutions  for  (5) »  say 
,  *Numbers  in  brackets  refer  to  references  listed  at  the  end 


>f  the  paper 


I 
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’/l  ’^2  f2 


taken  so  that  at  t  =  0  these  have  the  values 


1,0, 0,1,  respectively o  Then  every  set  of  solutions  Is  a  linear 
combination  of  these  two.,  Now  the  fundamental  equation  (6)  Is 
linear  and  homogenepus  In  the  pair  so  that  a  multiplicative 

factor  will  drop  out  from  the  formulas  for  6x,6y<,  Hence  we 
may  express  all  solutions  of  Interest  in  terms  of  a  single  param¬ 
eter  a 


(10) 


In  so  far  as  the  use  of  this  formula  is  concernedo  An  extremal 

which  starts  at  the  origin  at  t  =  0  is  characterized  by  the 

value  of  a  except  for  the  terminal  time  T,  and  the  problem 

C  <3  *^  Ofl^t  imuCrvO 

of  determining  “bourse  reduces  to  that  of  finding  a, To 

It  Is  shown  in  the  appendix  that  for  extremals  a  differential 
change  in  a  leads  to  a  differential  change  6x(T) ,6y(T)  in 
x,y  at'-  T 


where 


(12)  J 


A^iy^2  t=T 


If  both  a,T  are  changed  the  resulting  differential  changes 
' in  the  terminal  values  are 


(13) 


Ax(T)  =  x(T)6T  -  ,^(T)6a 
Ay(T)  =  y(T)6T  +  jA(T)6a 
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4.  Computational  routine « 

Let  us  guess  initially  values  for  a,T,  say  a^jTQ. 

Then  Integrate  simultaneously  the  original  equations,  a  funda¬ 
mental  set  solutions  to  the  adjoint  system  (6), 

with  determined  then  by  (10)  and  p  determined  by  (8); 

calculate  also  the  integral  Jo 

A  curve  thus  determined  will  not  go  to  the  desired  poln^ 

X,Y  generally,  but  to,  say,  ,Y^  ,  where  i  is  the  itera¬ 
tion  index.  In  equation  (13)  ,  set 
,  Ax.  =  X  -  X. 

(1^) 

lAy^^  =  Y  -  Y^ 

and  determine  6a^,6T^,  =  T^+bT^^. 

Continue  this  until  some  convergence  criterion  is  satis¬ 
fied,  say, 

(X  -  X^)^  +  (Y  -  Y^)^  <  £  , 

where  £  is  a  prescribed  number. 

Of  course,  if  a,T  are  guessed  too  far  from  the  correct  values 
it  may  not  converge;  this  has  not  been  a  problem  so  far. 

5.  General  end  conditions. 

In  the  most  general  case  we  may  have  a  function  g(x,y,t)  to 
be  minimized,  or  maximized  at  the  terminal  point,  subject  to  N 
constraints  of  the  form 

hj^(x,y,t)  =  0,  n  =  1,‘>'=°N. 

There  may  be  none,  one,  or  two  of  these  (N  =  0,1, or  2).  In  this 
case  the  cond^jtJ.ons  on  the  end  values ,  called  the  transversal  con¬ 
ditions,  which  must  be  satisfied  for  a  stationary  value  of  g  are 


]  (Nx3) 


-?- 


first  (see  Bliss  [(S] ,  pZ03) 

/dh^/dx  dh^/dy  dh^/dt' 

(15)  rank  (  5g/Sx  dg/dy  dg/Dt  >=  N+1 

I  >  ^  -^x 

and  second,  the  rank  of  the  matrix  obtained  by  omitting  either  of 
the  last  two  rows  must  also  be  N+lo 

The  first  N  rows  of  the  above  matrix  are  the  coefficients 
of  Ax,Ay>  6T  in  the  next  row  are  from  Ag,  and  the  last 

row  is  the  corresponding  set  of  coefficients  in  the  differential 
formula  (6)  rewritten  in  the  form 

TT- 

(16)  i'X^x  +/^Ay  -(^x-yyy)6T]^^^  =  A»Vp6p  dto 

In  the  problem  studied  first,  this  reduces  to  the  condition 
()ix  +yuy)^  0,  and  is  always  satlsfiedo 

One  way  to  determine  the  path  is  to  guess  a,T  as  before. 
The  correctional  routines  may  be  obtained  from 
Ah^  =  Sh^/dxAx  +  dh^/dyAy  +  dh^/6t6T, 

and  it  may  be  necessary  to  use 

c  a)v(T)  =  A(T)6T  +  (dX/da)^&a 


i 


Ayi/(T)  =  y(^(T)6T  +  (^/da)^6ao 
The  differentials  are  chosen  to  drive  residual  errors  to  zero. 

Example.  Consider  the  problem  of  getting  to  any  point 
where  x  assumes  a  specified  value  X,  in  minimum  time,  with 
no  constraint  on  y.  The  matrix  of  (15)  is  then 
/I  0  0 

0  0  1 


y  -()vxy/y)J 
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Hence  yu(T)  =  0.  This  yields  one  relation 

=  (-^^sin  a  a)6a  +  yu(T)6T  =  yuiT) 

for  6a, 6T;  the  equation  for  Ax  in  (13) 
is  the  other  equation  neededo 

An  alternate  procedure  is  to  continue  the  computation 
until  X  =  X;  this  determines  T^o  Then  obtained 

by  setting  ^(T)  ^^cos  "*’/^2  i®^^  ^i+1  “  seems 

that  the  only  ingenuity  required  is  in  guessing  starting  values 
and  in  setting  up  correction  routines  from  (15)  »  typical  comput-» 
ing  problems.  As  this  example  shows,  the  correction  routines  are 
not  unique  and  some  may  converge  better  than  others. 

6.  Comments. 

In  a  simple  example  requiring  about  100  time  steps,  each 
Iteration  took  about  one  second  and  a  path  was  obtained  in  from 
four  to  twenty  seconds  on  the  CDC  l604,  depending  on  the  function 
V  and  the  initial  guesses?  this  could  well  be  decreased  if 
desired.  On  the  other  hand,  the  use  of  empirical  data  and  the 
attendant  calculations  will  undoubted  increase  the  time. 

There  are  two  other  possible  methods  of  solution  of  this 
problem.  If  the  function  v  is  independent  of  t,  then  the 
order  of  the  system  may  be  reduced  by  eliminating  the  time  t. 

The  resulting  Euler  equation  is  of  order  two  and  may  be  solved 
by  a  relaxation  procedure  (Haltinerj,  Hamilton^  'Arnason  [3]),.  A 
limited  comparison  suggests  that  the  differential  correction 
procedure  converges  more  rapidly  if  the  number  of 

time  steps  is  large.  No  way  is  seen  to  extend  the  relaxation 
method  to  differential  systems  of  higher  order.  The  only  other 


-9- 


method  goes  by  various  names:  steepest  descent  or  ascent,  and 
gradient.  It  is  based  on  a  corrective  routine  similar  to  the 
constructive  proof  of  the  fundamental  lemma  of  the  calculus  of 
variations  (Courant  [l]  p  200) »  Discussions  suggest  that  all 
have  uncertainties  associated  with  convergence. 

Corner  condition.  If  there  are  any  corners  (points  of  the 
curve  where  v  is  discontinuous)  then  the  differential  formulas 

ti+ 

(6)  and  (15)  need  a  term  added  -[vCXcos  p  +yusin  p)]  6tj^ 

A 

for  each  point.  For  the  values  of  A,yu  associated  with  an  ex¬ 
tremal,  this  vanishes  (but  does  not  vanish  generally  if 
yu  =y(^i  f  etc.).  A  program  is  needed  to  check  for  this  condition, 
if  it  may  occur,  since  the  values  of  p(t^-)5  p(t^+)  are  not 
close  to  each  other.  Implicit  in  the  derivation  of  the  above 
correction  is  that  the  transients  Introduced  by  "coming  about" 
are  negligible. 

Many  important  conditions  which  must  be  checked  to  ensure 
that  the  resulting  curve  affords  a  minimum  must  be  neglected  in 
a  short  paper.  These  are  covered  completely , though  tersely, 
by  Bliss  ([6]  chapters  7,8).  Many  of  these  have  no  significance 
until  a  path  is  obtained. 

As  applied  to  shin  routing: ^  these  have  also  been  discussed 
by  de  Jong  [4P  ^contains  a  discussion  of  several  interest¬ 
ing  cases  of  particular  interest  in  air  navigation.  The  direc¬ 
tion  perpendicular  to  /\  defines  the  curves  of  constant  time 
associated  with  the  extremals  from  the  origin.  An  alternate 
method  of  determining  minimum  courses  is  to  determine  these 
isochrones,  as  done  by  Hanssen  and  James  [5] 5,  except  to  use 
extremals  and  the  transversal  conditions.  This  may  be  done 
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on  the  computer,  as  outlined  in  some  preliminary  reports 
[lO]  and  sections  8  and  10  of  this  paper. 

This  method  of  solution  is  an  application  of  the  methods 
which  Bliss  [8]  introduced  for  calculating  differentials  in 
ballistics  (see  also  Bliss  [9]>  Ch.  V  for  summary,  and  p.  125 
for  other  references). 
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Part  II  OPTIMUM  SHIP  ROUTING 
Introduction 

Various  other  aspects  of  optimum  ship  routing  are  taken  up  here.  The 
next  sections  treat  the  case  where  the  velocity  field  does  not  change  signifi¬ 
cantly  with  time,  and  may  be  considered  constant  over  short  periods.  The 
differential  equations  are  then  of  lower  order  and  the  problem  is  simpler. 
Two  methods  of  solution  are  discussed.  Each  makes  use  of  extremals, 
which  are  the  minimum- time  curves,  or  brachis trochrones ,  at  least  over 
short  distances.  These  may  be  used  together  with  a  relation  known  as 
'*trans ve rsality”  to  determine  equal-time  curves,  or  isochrones.  These  de¬ 
fine  the  boundary  to  the  points  which  the  ship  can  reach  at  any  particular 
time:  it  can  generally  reach  any  point  up  to  and  on  these  equal-time  curves, 
but  no  point  beyond  them  at  the  corresponding  time;  also  points  on  the  iso¬ 
chrones  can  be  attained  only  by  extremals.  The  shortest-time  route  maybe 
determined  from  these,  A  second  method  of  solution  is  given  for  finding 
the  shorte st- time  route  directly  without  using  the  isochrones. 

Next  the  problem  of  determining  the  equal-time  curves  is  taken  up  for 
the  case  where  the  velocity  field  varies  rapidly  or  significantly  with  time. 
Then  some  routines  are  given  for  effecting  rendezvous  between  two  ships, 

(a)  when  one  ship  is  following  a  known  course,  and  (b)  when  the  two  ships 
cooperate.  Finally  a  brief  discussion  is  given  of  problems  which  may  be 
enccuntevel  in  the  computation. 


-13-  Pa  #32 

/  Constant  Velocity  Field,  Euler  Equations,  Transversals. 

In  this  section  It  will  be  assumed  that  the  speed  v  Is 
a  function  of  position  and  heading,  but  not  time,  v  =  v(x,y,y' 
where  x,y  are  position  coordinates,  and  y'  =  dy/dx.  For 
simplicity,  we  will  generally  choose  the  coordinate  set  so  that 
the  Initial  point  Is  the  origin  and  the  final  point  (X,0) 

Is  on  the  x  axis. 

The  Euler  equation,  which  Is  a  necessary  condition  for 
minimum  time  will  be  derived  now.  If  x,y  are  cartesian  coord 
Inates,  the  time  required  to  go  a  short  distance  along  a  curve 
Is  approximately 

(17)  t  =  As/v  =  /l+y'^  x/v  , 

and  the  time  to  go  along  any  selected  curve  Is  then 


(18) 


T  = 


X 


f’(x,y  ,y'  )dx; 


1  2 

this  defines  f  =  -Nll+y'  /v.  If  x,y  are  not  cartesian 
coordinates ,  there  are  corresponding  relations ,  depending  on 
the  metric  of  the  coordinate  set. 

Let  us  consider  two  neighboring  paths.  On  the  second, 
y-^y  +  6y 


(19) 


y'  ■  >  y'  +  6y'  , 


where  (6y)'  =  6(y' ) ,  and  means  "Is  replaced  by".  The 

f  I 

Integral  |6y'l  dx  Is  "small".  The  difference  In  the  time 

^0 

required  to  follow  the  two  courses  may  be  approximated  then 

/  X 


(20) 


6T  = 


(fy&y  +  fyi&yMdx  +  f  (x,y ,y' )^^J^6X  ; 


subscripts  In  the  Integrand  indicate  corresponding  partial 
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derivatives.  Let  us  integrate  this  by  parts  to  eliminate  &y' 
from  the  integrand 

("X  X 

(f„  -  df  , /dx)6y  dx  +  f(x,y,y*)Y&X  +Lf^,6y]  . 

0  y  y  ^  y  0 

If  the  end  points  are  fixed,  then  6X  =  6y(0)  =  6y(X)  =0.  If 
T  is  a  minimum,  because  of  the  choice  of  path,  then  6T  must 
vanish  for  all  allowed  functions  6y.  Since  6y  is  arbitrary 
except  for  being  "small"  and  vanishing  at  the  endpoints  of  the 
curve,  it  follows  that 


-  fy  =  0- 


This  is  the  well-known  Euler  equation  and  is  a  necessary 
condition  (assuming  that  f  has  continuous  second  derivatives 

We  may 

rewrite  this  as 


and  f'yiyi  5^  0)j  see  Bliss  [6],  p.  Ilf  for  a  proof 

»/  J 


(23) 


y"  = 


fy,y,^^y  "  ^'^y'y  ^y'x-^* 


which  is  of  the  form  y"  =  P(x,y,y').  It  may  also  be  rewritten 
as  a  pair  of  first-order  equations,  by  introducing  a  new  variable 


(24) 


y'  =  Z 


z'  =  F(x,y,z) , 

which  is  convenient  for  computations.  It  should  be  observed 
that  since  (23)  is  of  second  order,  there  are  two  constants  of 
integration  and  an  extremal  is  determined  if  the  values  of  y,y* 
are  given  for  some  value  of  x;  since  y(0)  =  O.ln  our  problem, 
there  is  a  single  constant  of  integration  to  be  determined. 

The  curves  defined  by  (22)  or  (23)  are  called  extremals . 
Transversals.  The  equations  for  curves  which  correspond 
to  equal  times  along  different  extremals  of  a  family,  such  as 
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the  family  passing  through  a  given  point,  will  be  derived  now. 
If  the  endpoint  (X,Y)  is  not  fixed,  but  may  vary,  then 


the  terms  in  equation  (21)  associated  with  it  do  not  necessarily 
vanish.  Since  we  are  interested  in  shortest-time  routes,  we 
will  consider  only  extremals,  so  that  the  integrand  and  the  inte¬ 
gral  vanish.  Now  the  change  in  the  final  value  of  y  is 


(25) 


6Y  =  6y(X)  +  y' (X)6X. 


i6y(X) 


6Y 


(see  sketch),  y'  being  the  value 
on  the  extremal.  Equation  (21) 
becomes 

(26)  6T  =  [f6X  -t(£Y  -  y'6X)fy,]j^ 


J  ly'6xy 

6X 

Relations  among 
variations  of 
end  values. 


Now  if  we  have  a  one-parameter  family  of  extremals  such 


as  those  coming  from  a  given  point,  then  the  various  points, 
one  on  each  extremal ,  which  correspond  to  the  same  value  of  T 
will  be  obtained  by  setting  6T  =  0.  That  is,  from  (25)  and 


(26), 

(27)  it  -  [y’  -  t/t  ,]^bx  =  0; 

this  is  the  defining  relation  for  the  transversal  direction 
6Y/6X.  A  curve  S  which  cuts  each  of  a  family  ^e]-  of  extre¬ 
mals  transversally  is  transversal  to  the  family,  or  is  a 
transversal  of  the  family.  Its  equation  is 


(28)  dY/dX  =  y'  -  f/fy,  ; 

dY/dX  is  the  slope  of  S  and  y'  is  the  value  on  the  cor- 

c. 

responding  extremal. 

And  if  we  are  given  a  smooth  curve  S,  then  there  is  a 
family  of  extremals  |E|  which  satisfy  (28).  The  family 
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so  defined  represents  the  curves  of  shortest  time  from  the 
transversal  S,  at  least  if  the  distances  are  not  too  great. 


Comment.  Sometimes  the  transversals  are  perpendicular 
to  the  extremals.  It  may  be  easily  verified  that  this  is  the 
case  for  an  arbitrary  family  of  extremals  if  and  only  if  f 


in  the  ship-routing  problem 


has  the  form 


this  is  the  case  if  v  is  independent  of  the  heading. 

8.  Numerical  routine. 

A  method  of  determining  a  minimum-time  ship  course  is 
now  given.  It  is  similar  to  the  method  described  by  Hanssen 
and  James  [5]  (pp259»260)  except  that  it  makes  use  of  the 
results  Just  derived,  Involving  extremals  and  transversals. 

Let  us  take  the  set  of  extremals  which  emanate  from  the 
starting  point.  Calculate  a  family  of  these  with  the  initial 
heading  angle  as  a  parameter,  say  taking  values  of  ban”^y'(0) 
one  degree  or  ten  degrees  apart,  using  the  Euler  equation. 
Qontinue  these  for  say  six  or  twenty-four  hours.  The  value  of 
t  must  be  obtained  by  integration  and  the  terminal  point  will 
involve  interpolation.  At  the  endpoints  determine  also,  from 
(28)  the  transveral  direction. 

We  now  have  a  set  of  points  on  the  equal-time  curve  and 
the  corresponding  direction.  Fit  a  curve,  .  using  this 

data.  It  is  not  the  most  common  type  of  data  since  not  only 
is  the  point  given,  but  also  the  slope.  Name  this  curve  , 
or  case  may  be.  It  represents  the  maximum  distance 
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the  ship  can  be  at  that  time,  based  on  the  Initial  velocity 
field. 

Now  revise  or  update  the  velocity  field  to  correspond  to 
the  values  forecast  for  six  or  twenty-fours  ahead.  Then  take 
the  one-parameter  family  of  extremals  which  cut  ,  or 

transversally.  The  slopes  of  each  extremal  is  given  by  equation 

(28) ,  with  revised  values  for  v,  for  the  continuing  extremals. 
As  many  extremals  as  desired  can  be  drawn  out  from  S^.  With 
each  extremal  that  is  continued  from  associate  the  corres¬ 
ponding  value  of  y'(0). 

Continue  these  out  for  another  six  or  twenty-four  hours 
and  generate  another  transversal  ^48*  Update  the 

field  and  continue  until  a  transversal  hits  the  desired  terminal 
point.  Interpolate  to  get  the  initial  heading  and  the  route. 

This  method  has  the  desirable  feature  that  the  results 
may  be  easily  Interpreted,  particularly  by  those  already  cal¬ 
culating  minimum  time  routes  by  present  methods.  On  the  other 
hand  it  requires  unnecessary  computing  to  determing  one  route. 

9.  Alternate  numerical  routine. 

It  is  felt  that  the  following  method  will  generally 
determine  a  route  more  quickly. 

Let  us  consider  equations  (24).  The  equations  for  the 
variations  of  an  extremal  are 
6y'  =  6z 

(29) 

6z'  =  Fy6y  +  F^6z, 

Let  us  guess  an  initial  value  for  y' ,  Then  compute  the 
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corresponding  extremal  ,  which  will  probably  not  go  through 

the  desired  endpoint  (X,0).  Simultaneously,  compute  a  solution 
to  (29) »  with  6y(0)  =  0  and  any  convenient  value  of  6z(0)  , 
since  the  system  is  homogeneous  and  linear  in  6y,6z<,  The 
value  of  y(X)  will  be  in  error,  say  it  is  Treat  6y(X) 

as  a  difference;  set  6y(X)  =  -  and  solve  for  6z  =  6y'(0). 
This  gives  a  corrected  value  for  y'(0).  Continue  this  until 
the  terminal  value  differs  from  the  desired  value  of  y  by 
less  than  some  preassigned  number,  the  convergence  criterion. 

The  values  for  the  velocity  field  could  be  updated  regularly 
every  six  hours  as  before,  if  desired.  A  limited  number  of 
computations  suggests  that  this  will  give  a  solution  quickly. 

Comment.  There  is  a  condition  for  an  extremal  to  yeld 
a  minimum  value  for  the  integral,  called  the  non-tangency  con¬ 
dition.  If  at  any  point  &y  as  obtained  from  (29)  above  is 
zero,  the  corresponding  point  is  on  the  envelope  of  the  family 
of  extremals  and  does  not  furnish  a  minimum  value  to  the  integral. 
In  this  case  a  warning  should  be  given  to  the  operator,  or 
a  corresponding  subroutine  generated.  This  will  probably 
never  occur  for  short  courses. 

10.  Isochrones  for  time-varying  field. 

A  method  is  given  here  for  constructing  the  Isochrones 
when  the  velocity  field  varies  with  time,  and  for  determining 
the  minimum  time  course  by  the  method  of  section  8. 

Let  us  consider  various  extremals  starting  from  the  origin; 
each  is  determined  by  a  value  of  a  in  equation  (10)  ,  and 
by  equations  (1),(5)>(8).  At  every  point  of  each  extremal  E 


Interrelations  between 
extremal  E,  velocity  v, 
solution  A  to  adjoint, 
and  the  Isochrone  So 


the  relation  between  the  velocity  v,  the  solution  A  to  the  ad¬ 
joint,  the  extremal  E  and  the  Isochrone  S  is  as  indicated  in 
the  figure o  Let  r*  be  the  curve  defined  by  the  allowed  velocity 
V  at  P,  with  p  as  a  parameter. 

Then  the  heading  p  is  chosen  so 
that  V  has  a  maximum  projection 
onto  A  ,  see  section  2,  The 
curve  S  perpendicular  to  A  at 
P  is  generated  by  small  changes 
in  the  parameter  a;  it  is  a  trans¬ 
versal  of  the  family.  The  impor¬ 
tant  property  of  S  is  that  all 

points  on  it  and  to  its  left  in  some  neighborhood  can  be  reached 
at  that  time,  no  points  to  the  left  can  be,  and  points  of  S 
are  attained  only  by  extremals.  That  is,  S  is  part  of  the 
boundary  of  a  closed  region  whose  points  are  exactly  the  points 
where  the  ship  may  be  at  that  time  t. 

We  then  get  points  on  a  curve  S  and  the  tangent  direction. 
These  define  the  isochrone  S(t),  A  relatively  small  number  of 
points  is  required  to  define  S  since  the  tangent  direction  is 
also  given;  this  tends  to  be  offset  by  the  fact  that  v  is  given 
from  empirical  data.  If  t  is  small  enough,  the  curve  S  is 
similar  to  a  circle  or  an  arc  thereof.  If  the  isochrone  S  is 
given,  it  must  be  emphasized  that  the  normal  to  it  determines 
A  ,  and  that  the  extremals  are  not  generally  perpendicular  to  S, 

The  first  isochrone  S(t)  which  touches  the  specified  terminal 
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point  obviously  yields  the  minimum-time  route.  The  value  of 
T  and  the  parameter  a  must  be  determined  by  Interpolation, 

If  the  total  time  exceeds  the  time  that  forecasting  can  be 
done  reliably,  then  the  extremals  may  be  replaced  by  terminal 
segments  of  great  circle  routes,  the  extremals  for  a  uniform 
velocity  field,  as  done  by  Hanssen  and  James  [5]  (p  26I)  ,  or 
statistical  data  may  be  used  as  the  basis  for  determining  v. 
Whatever  data  is  "most  reliable"  should  be  used;  the  results 
can  be  no  more  reliable  than  the  data, 

11.  Rendezvous  between  ships. 

Two  ships  will  be  said  to  rendezvous  if  at  some  time 
their  positions  coincide.  It  is  assumed  that  the  time  required 
for  the  terminal  maneuvering  is  negligible  compared  to  the 
total  time  spent. 

In  the  first  example,  let  the  position  of  the  second 
ship  be  denoted  by  x*,y*,  these  being  known  functions  of 
time.  As  before  (see  section  4),  we  guess  values  for  the 
parameters  of  the  extremal,  a,T.  In  the  correction  routine 
for  these  we  must  allow  for  the  distance  the  second  ship 
will  travel  if  we  change  T  by  an  amount  6T.  In  place  of 
equations  (13),(l4)  ,  use 

X*  -  X  =  (5;  -  x*)6T  - 

(29)  .  1 

y*  “  y  =  (y  -  y*)6T  +')\J5a, 

for  corrections  6a, 6T  to  the  estimates  of  a,T,  all  quan¬ 
tities  in  (29)  being  evaluated  at  t  =  T  after  the  first 
estimate  of  the  trajectory  has  been  calculated,  J  was  defined 
by  equation  (12). 
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The  corrected  values  for  a,T  are  used  and  another 
course  computed.  This  Is  continued  until  the  distance  be¬ 
tween  the  two  is  less  than  some  preassigned  number, at  the 
end  of  a  computation. 

Cooperating  ships.  The  problem  is  more  interesting 
if  the  second  ship  cooperates.  Let  quantities  associated 
with  the  second  ship  be  denoted  by  an  asterisk  (*).  For 
rendezvous,  there  must  be  a  time  T  such  that  x(T)  =  x*(T), 
y(T)  =  y*(T).  There  is  a  further  condition  that 

(30)  A(T)  II  A*(T)  ; 

the  vectors  defined  by  the  adjoint  systems  must  be  parallel 
at  time  T.  That  is,  if  H  =  then 

(31)  H  =  0. 

We  will  need  also  the  differential  change  in  in 

our  computational  routine.  If  we  change  a,T  by  small  amounts, 
we  get  a  differential  change  in  the  terminal  values 

f 

aA  =  [-Aj^(T)sln  a  +7\2(T)cos  a]6a  +  A(T)6T 

(32) 

=  [•^j^(T)sln  a  aj6a  +  y>'(T)6T 

and  two  more ,  exactly  like  these,  with  starred  terms,  asso¬ 
ciated  with  the  second  ship. 


Computational  routine.  Guess  a,a*,T  and  calculate  a 
first  approximation  to  the  courses  for  the  two  ships.  Let 
the  values  at  the  end  of  the  I'th  iteration  be  x^^  ,yj^  ,x*  ,y* . 
Generally  x^^  x* ,  or  y^^  y* ,  or  A  is  not  parallel  to  A" 
at  t  =  T,  Use  the  differential  formulas  as  differences  to 
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drive  the  errors  or  residuals  to  zero: 


(33) 


=  (Xj^  -  x^)6T  -  J^ba  +  J^*6a* 
y*  -  =  (yj^  "  y*)6T  +  j'Xba  -  J*'A*6a^ 


-(y* 

+  [(-A]^sln  a  a^* 

-[(-P\*sin  a*  +;>^2COS  a*J^ 


-  pM  -y^X*)6T 

-(70'lSin  a  yj^cos  a)A*]&a 

-(y^*sln  a*  yJ^Qos  a*)'}\  ]6a*. 


If  this  converges,  it  yields  a  stationary  value  of  T.  The 
computation  would  be  continued  until  some  convergence  criterion 
is  satisfied,  say  until  (x^-x*)^  +  (y^^-y*)^  +  <  £  ,  where 

c  is  a  preassigned  number. 

An  alternate  procedure  is  the  following.  Suppose  the  co- 
orlnates  are  chosen  so  that  the  first  ship  is  initially  at  (0,0) 
and  the  second  at  (X,0),  with  X  >  0.  Guess  a, a*  as  before 
and  compute  the  two  extremals ,  stopping  when  at  some  time  t  = 
T^,  X  =  x*.  Equations  (33)  again  constitute  a  set  of  three 
equations  for  three  unknowns,  though  only  two  of  these,  6a, 6a*, 
are  needed  to  start  the  next  iteration. 

An  alternate  method  would  be  to  determine  the  Isochrones 
S(t)  ,S*(t)  for  each  ship,  by  the  method  described  in  section 
10,  continuing  until  they  touch.  Since  (or  if)  they  are  smooth 
curves ,  S  and  S*  are  tangent  to  one  another  at  the  first 
time  of  contact.  This  is  equivalent  to  condition  (31) »  that 
AHA*,  since  A  and  S  are  alvjays  perpendicular. 
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12.  Comments 

As  remarked  initially  the  first  purpose  of  this  paper  has  been 
to  give  a  method  for  determining  optimum  ship  routes..  It  may  have 
given  however  an  impression  that  the  problem  is  simpler  than  it 
actually  is.  In  this  section  a  brief  discussion  is  given  of  some 
situations  which  may  arise,  which  are  not  treated  earlier. 

It  should  be  remarked  that  the  routine  given  converges  to 
a  stationary  value  of  time,  which  may  not  be  a  minimum.  This  is 
particularly  the  case  for  long  routes,  where  a  course  not  near 
that  chosen  may  yield  a  smaller  value  of  time.  There  is  some  con¬ 
dition  which  is  an  extension  of  the  envelope  condition  in  elemen¬ 
tary  calculus  of  variations,  and  the  well-known  conditions  of 
Weierstrass  and  Legendre  (see  Bliss  [6],  chapter  8).  This  does 
not  lessen  the  value  of  this  paper,  since  these  conditions  can 
be  checked  only  after  a  course  has  been  determined. 

The  envelope  condition  in  the  simplest  case  states  that  if 
an  extremal  furnishes  an  extreme  value  to  an  integral,  then  it 
cannot  contain  a  focal  point  or  a  point  of  the  envelope  of  the 
family  of  extremals  from  the  initial  point.  For  the  simplest 
problems  this  can  be  checked  by  checking  to  see  if  6y  is  zero 
at  any  point  of  the  course,  but  no  simple  check  is  known  generally 
to  the  author. 

If  the  course  is  determined  by  the  method  of  isochrones,  then 
it  is  suspected  that  the  Isochrone  will  then  develope  an  interior 
corner.  Unfortunately  the  method  of  fitting  the  Isochrones  will 
probably  have  an  implicit  assumption  of  smoothness  which  may  well 
obscure  the  existence  of  the  corner.  The  difficulty  is  aggravated 
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by  the  fact  that  the  data  is  given  empirically  and  we  must  try 
to  smooth  out  or  minimize  the  random  errors  by  a  curve-fitting 
or  smoothing  routine. 

For  the  initial  part  of  the  course  we  will  surely  be  able 
to  predict  accurately.  For  the  later  part,  we  may  be  required  to 
use  statistical  data,  or  in  the  absence  of  that,  a  terminal  sec¬ 
tion  consisting  of  part  of  a  great  circle.  If  the  times  involved 
are  large,  so  that  better  data  becomes  available,  the  course 
should  be  continually  redetermined,  using  the  current  position 
as  tile  starting  point.  This  is  particularly  important  on  long 
routes,  and  when  storms  Introduce  large  changes  in  the  speed 
which  cannot  be  estimated  accurately  beforehand. 

It  seems  that  time  may  well  be  too  simple  a  cost  function. 
Cargo  damage,  danger,  passenger  discomfort,  etc.,  increase  with 
wave  height.  The  cost  function  must  be  made  up  by  someone 
familiar  with  the  details  of  shipping.  In  the  best  routing 
procedures,  several  activities  must  be  coordinated.  (1)  First, 
the  ^meteorologist  or  oceanographer  must  collect  the  data  and 
predict  sea  state,  and  weather,  too,  if  it  is  significant. 

(2)  A  cost  engineer  familiar  with  the  ship  and  its  cargo  must 
make  up  a  realistic  cost  function  in  terms  of  time,  danger, 
etc.  (3)  The  mathematical  framework  is  given  here.  (4)  Finally, 
there  is  no  little  problem  of  programming  properly.  Some  study 
needd  to  be  made  of  methods  of  smoothing  the  data,  so  that  the 
derivatives,  which  are  notoriously  poor  when  the  data  is  poor, 
are  satisfactorily  smooth.  The  success  of  the  method  depends 
on  the  Integrated  efforts  of  all  of  these. 
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Finally,  let  us  consider  a  problem  where  the  above  routine 


breaks  down.  Let  us  consider  a  sailboat  In  a  narrow  channel 
where  hills  on  either  side  slow  down  the  wind.  Let  us  consider 


p  p 

speed  In  the  form  v  =  v^(l  -  y  /w  )(1  -  2  cos  p) ,  where 


and  w  are  constants.  It  Is  easily  verified  that  the  maximum 
speed  upwind  (In  the  direction  of  the  positive  x  axis  Is  v^/4. 
However,  to  effect  this,  the  boat  must  come  about  an  Infinite 
number  of  times;  this  Is  known  In  control  theory  as  the "chattering" 
solution.  The  above-given  routine  is  too  simple  for  this  problem: 
the  time  spent  coming  about  and  the  ensuing  transients  are  not 
negligible.  It  seems  to  the  author  that  some  feel  for  the 
mechanics  of  such  problems  is  necessary  in  the  programmer,  that 
one  must  suspect  beforehand  that  such  events  are  likely.  These 
seem  to  be  of  little  importance  with  powered  ships  and  were  not 
investigated  in  any  detail. 
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Appendix.  In  this  section  the  differential  formulas 
of  section  3  are  derived.  Let  Ai,A2  denote  as  before  the  so- 
lutions  to  the  adjoint  system  which  have  initial  values  i,j, 
resp.  Let  a  =  a^cos  a  +  AgSin  a,  as  before  and  a  =  |a|. 

-4  -4 

Then  =  Aj^cos^a  +  Ag^sin^a  +  2  A^'Agsin  a  cos  a.  Let 
X  =  Vcos  Pj  y  =  V  sin  p,  and  denote  the  "velocity"  by  V; 
lower  case  \^v  •  was  used  earlier.  Lower  case  symbols  will  be 

—4  — • 

used  to  denote  unit  vectors  in  this  section,  V  =  Vv  ,  a  = 

AW,  Then 

7  =  v^v  +  Vn, 

p  p  » 

^4  »4 

where  n  =  k^v.  If  u  is  the  unit  vector  parallel  to  V^, 
then 

u  =  k  X  a/a 


and  - 

V  =  /(V  2+V2)  u 

-4  -4 

on  every  extremal,  since  V^^a  s  0. 

Now  let  us  consider  at  the  same  time  and  place  two  ex¬ 
tremals  differihg  by  small  amounts  6p,6ao  Then 


''pp-''  5P  +  %.Ag6a  =  0. 

Now  also 


So 

(A.2) 


IVpIpU  IVpIkxu 

iVpIpk!^-  !Vyiw> 

\ 

=  -/(V2+Vp2)  A. 


^4 

^  II  ‘AW 


Now  consider  also 


^The  terms  to  the  right  of  the  parallel  bars  ||  in  the 
right-hand  margin  will  be  used  to  suggest  the  operations  by 
which  the  next  equation  is  obtained. 


»Ao2- 


V  »«A  =  0  =  V  »(aiCOs  a  +  Aosin  a) 
P  P  ^ 


Let  also 

P*  = 

(Vp-Al)®+  (Vp-Xj)2  j 

then 

sin  a 

=  -Vp*Ai/p  ,  oos  a  =  Vp-Aj/p 

and 

(A.3) 

V  .A 

P  a 

=  Vp»(-Aisin  a  +  Agcos  a) 

=  Po 


Hence  from  (Ad)  ,(A,2)  ,(A„3) 

6p  =  p6a/(Ay[V2+Vp2]), 

and  the  differential  formula  (6)  derived  earlier  becomes 
when  A  =  Ai 

»T 

(Xi6x  +  |ii6y)m  =  -  (p^sin  a/A)dt  da 

1  Jq 

for  extremals.  But 


p2  =  [y(V2+Vp2)^-Aj^  +  L/(V2+Vp2)2^»A2]* 

=  ^  {[(Ajcos  a  +  Agsin  a)  "Aixjtjz 

a2  '■ 

+  ['(a^cos  a  +  Agsin  a)«A2xIc]2| 


Hence 


and 


V2+V  2 


Xi  Xg 
dl  U2 


I 

(Xi6x  +  ^^6y)^  =  -J 


(X26X  +  M26y)^  =  J 


-J 

pT  y2^y  2 

Xi  X2 

'0  A^ 

Ml  M2 

=  J 

pT  y2+y^2 

Xi  X2 

'0  A^ 

Ml  M2 

dt  sin  a  6a  , 


l|-d2 


dt  cos  a  6a  , 


II  Ml 


2 
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Hence 


6x(T)  =  -  J  w(T)  6a 
6y(T)  =  J  X(T)  6a, 

the  values  of  X,u  being  those  associated  with  the  maximizing 
condition  and 


J  = 


^2 

Ui  Ug  1 


V2+V  2 

Xi,  Xg 

‘'O 

dt 


This  derivation  seems  long  and  complex  but  no  simplification 

has  been  found,  and  the  corresponding  three-dimensional  formula 

has  so  far  defied  derivation.  ^ 

Frank  D,  Faulkner  6  June  1962 
Monterey,  Calif. 
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